Spreading of wave packets in disordered systems with tunable nonlinearity 
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We study the spreading of single-site excitations in one-dimensional disordered Klein-Gordon 
chains with tunable nonlinearity for different values of a. We perform extensive numerical 

simulations where wave packets are evolved a) without and, b) with dephasing in normal mode 
space. Subdiffusive spreading is observed with the second moment of wave packets growing as t". 
The dependence of the numerically computed exponent q on cr is in very good agreement with our 
theoretical predictions both for the evolution of the wave packet with and without dephasing (for 
cr > 2 in the latter case). We discuss evidence of the existence of a regime of strong chaos, and 
observe destruction of Anderson localization in the packet tails for small values of a. 
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I. INTRODUCTION 

The presence of uncorrelated spatial disorder in one- 
dimensional linear wave equations results in the local- 
ization of their normal modes (NMs). This is the well- 
known phenomenon of Anderson localization ^Ij and has 
been experimentally observed in a variety of systems, in- 
cluding as examples light 0, Q and matter [J ^ waves. 

Understanding the effect of nonlinearity on the local- 
ization properties of wave packets in disordered systems 
is a challenging task, which, has attracted the attention 
of many researchers in recent years. Several numerical 
studies of wave packet propagation in different models 
showed that the second moment m2 of norm/energy dis- 
tributions grows subdiffusively in time following a power 
law of the form m2 |6l-[ll|. 

In [ll| the mechanisms of spreading and localiza- 
tion were studied for the disordered discrete nonlinear 
Schrodinger equation (DNLS) 



(1) 



and the quartic Klein-Gordon lattice (KG) of anharmonic 
oscillators with nearest neighbor coupling. The exponent 
a was numerically found to be close to a w 1/3 and a 
theoretical explanation of this value was provided. 

Applying the same theoretical argumentation to a gen- 
eralized DNLS (gDNLS) model 
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with a being a positive integer, the dependence of a on ct 
was predicted in Q . The validity of this estimation has 
not been analyzed in detail. Mulansky [12| presented 
numerical simulations of the gDNLS model for a few in- 
teger values of a. In [isj numerical simulations of the 
gDNLS model were performed for non integer values of a 
on rather short time scales, leaving the characteristics of 
the asymptotic (t — > oo) evolution of wave packets aside 
and open. 

The main scope of the present work is to verify the 
validity and the generality of the theoretical predictions 



presented in Q for one-dimensional disordered nonlin- 
ear chains as a function of a. In particular, we choose 
to perform numerical simulations for the generalized KG 
(gKG) model instead of the gDNLS system. This choice 
was done for two reasons. Firstly, it allows us to test 
whether the estimations obtained in hold irrespec- 
tively of the presence of a second integral of motion (the 
norm J^i IV'iP for the gDNLS model). The second rea- 
son is a practical one. From the comparative study of 
DNLS and KG models (which correspond to ct = 2) per- 
formed in 0, [ll| it was observed that the KG model 
requires less CPU time than the DNLS system in order 
to be integrated up to the same time with the same pre- 
cision. Since in our study we are mainly interested in 
the characteristics of the asymptotic dynamical behavior 
of wave packets, the gKG model was preferred, as it per- 
mits long integrations of large lattice sizes within feasible 
CPU times. 



II. THE GENERALIZED KLEIN-GORDON 
MODEL 



The Hamiltonian of the gKG model is 
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where I is the lattice site index, ui and pi are respec- 
tively the generalized coordinates and momenta, a de- 
fines the order of the nonlinearity, W denotes the disor- 
der strength and ei are chosen uniformly from the interval 
[i, |] . The case ct = 2 corresponds to the standard KG 
model, which exhibits a similar dynamical behavior with 
the DNLS model [ll|. The equations of motion are 
ill = —dH/dui and yield 
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Hamiltonian ([3]) is a conservative system and its total 
energy E' > is preserved and serves as a control param- 
eter of the nonlinearity strength. Equations become 
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linear by neglecting the nonlinear term ImjI'^u; or in the 
limit ^ ^ (i.e. \ui\ 0) where < \ui\. Then 

the ansatz ui = Ai exp{iujt) reduces them to the linear 
eigenvalue problem 

c.^A.= (e+A)A.-(i,)(A_,+A,+r). (5) 

The normalized eigenvectors Ai,^i {J2i^li = 1) ^^'^^ 
NMs of the system with the corresponding eigenvalues 
Ai, = uj'^. The width of the squared eigenfrequency 
spectrum is Ak = I + with ujI £ f + ■ 

The asymptotic spatial decay of an eigenvector is given 
by A^^i - e-'/«(^-) where ^(X^) < C(O) ~ is 
the localization length • The spatial extend of a NM 
can be characterized by the average localization volume 
V ~ •\/12to2, with 7712 being its second moment [l5l |. 
F « 1 for > 10 , and F « 3.6^(0) for W < A. 
The average spacing of squared eigenfrequencies of NMs 
within the range of a localization volume is d = /V ~ 
_^ 4iy)/360 for < 4 and d - A^f for > 10. 

The squared frequency shift of a single-site oscillator 
induced by the nonlinearity is 

f E \ ^ 2^ r (— ) 

'"'^ = 0F(a+2)r(4^) ' 

where Ei is the energy of the oscillator (see Appendix 
For o- 2 in Eq. © it follows 5i « {■iEi)/{2ii) 

In our study W = A and we follow the evolution of 
single site excitations for long times and for < cr < 4. 
Then we have A^ = 2, V" w 20 and d w 0.1. We excite 
site Iq by setting pi = ^/2E5i,ig, m; = for t = with 
ejo = 1. In our simulations we use symplectic integration 
schemes of the SABA family of integrators [ll|, [la| with 
some particularities [T7| . 

In our computations the number of lattice sites N and 
the integration time step r varied between N — 500 to 
N = 3000 and t = 0.2 to t = 0.05, in order to exclude 
finite size effects in the wave packet evolution and allow 
long integrations up to 10^ time units. In all our sim- 
ulations the relative energy error was kept smaller than 
10-3. 

Following the methodology of 0, [ll| we order the 
NMs in space by increasing value of their center-of- 
norm coordinate AT^ = consider normal- 

ized energy density distributions z^, = Ei,/ E^^ with 

El, = Al/2 + ujIAI/2, where A,j is the amplitude of 
the vth NM and the corresponding squared eigen- 
frequency. In our analysis we use the second moment 
TO2 = ~ v^Zi, {v = ^^vz^), which quantifies 

the wave packet's degree of spreading, the participation 
number P — l/X^jy-^^i which measures the number of 
the strongest excited modes in z^, and the compactness 
index Q = P^/m2, which quantifies the sparseness of a 
wave packet, since C decreases as the wave packet be- 
comes more sparse (see [ij] for more details). 



III. DIFFERENT DYNAMICAL REGIMES 

DiiTerent dynamical regimes have been discussed in re- 
cent publications 0-S El • In a recent paper one of us 
presented a coherent interpretation of collected numeri- 
cal data within a simple framework which uses two av- 
eraged parameters of an initial wave packet as essential 
control parameters for the dynamical evolution: the av- 
erage norm or energy density in a packet, and its typical 
size [l^. According to tjiese results, a wave packet can 
be selftrapped (see also i3,|lli]) in a regime of strong non- 
linearity. This happens when nonlinear frequency shifts 
are larger than the width of the spectrum of the linear 
equations. If the nonlinearity is weak enough to avoid 
selftrapping, the wave packet will spread either in an in- 
termediate regime of strong chaos, followed by an asymp- 
totic regime of weak chaos, or the strong chaos regime 
is skipped, and spreading starts in the regime of weak 
chaos and stays there. The outcome depends on the ra- 
tio of the nonlinear frequency shift of the wave packet 
after spreading into the full localization volume and the 
average spacing d. Here we will adapt these arguments 
to the study of single site excitations. Then the three 
possible regimes are: 

a^E"/'^ > Ak : selftrapping, (7) 
Oct ( — j > d : strong chaos, (8) 

Oct f — j < d : weak chaos. (9) 

Wave packets in the strong chaos regime spread faster 
than in the weak chaos case (see also Sect. lIVAl below). 

The location of the three different dynamical regimes 
in the parameter space of the system's energy E and the 
order a of the nonlinearity is shown in Fig. [1] For cr > 2 
the regime of strong chaos is absent (in the sense that if 
at all, it will coexist with selftrapping). Therefore, if not 
selftrapped, the wave packet is expected to spread in the 
asymptotic regime of weak chaos. Following Anderson's 
definition of localization 1], we measure the fraction Ey 
of the wave packet energy in a localization volume 1^ = 20 
around the initially excited site. For a localized state this 
fraction should asymptotically tend to a nonzero con- 
stant. We find that in the weak chaos regime (lower 
inset of Fig. [T]) the fraction continuously drops down in 
time, indicating complete delocalization. Contrary, in 
the strong nonlinear regime of selftrapping (upper inset 
of Fig. [ij the fraction appears to tend towards a nonzero 
constant. These behaviors are also clearly seen in Fig. [21 

For (7 the selftrapped regime is shifted to very 
large energies. At the same time for cr < 2 the regime 
of strong chaos is widening its window with decreasing 
CT. A wave packet, if not selftrapped, may then spread in 
the intermediate regime of strong chaos, and cross over 
to weak chaos at energy densities which are the smaller 
the smaller a is. Therefore the crossover to the asymp- 
totic regime of weak chaos may be pushed to very large 
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FIG. 1: The three different dynamical regimes for the gKG 
model, in the parameter space of the nonhnearity order ct 
and the energy i5 of a single nonlinear oscillator excitation. 
The points correspond to the particular numerical simulations 
presented in Sect. I1V"B] Insets: The fraction Ev of the wave 
packet's energy in a localization volume V = 20 around the 
initially excited site Zo versus time t in log- log plots for cr = 2.5 
and total energy E = 0.45 (weak chaos regime, lower inset), 
E = 2.0 (selftrapped regime, upper inset). The disorder real- 
ization is the same as in Fig. [S] 





FIG. 3: (Color online) Different spreading behaviors. m2 and 
P versus time t in log-log plots for different nonlinearity orders 
(7. Left plots: a = 1.5 and total energy E — 0.03,0.35,2.2 
[(b), blue; (g) green; (r) red]. Right plots: a = 2.5 and total 
energy E = 0.11,0.45,2.0 [(b) blue; (g) green; (r) red]. The 
disorder realization is the same for both values of a. Straight 
lines guide the eye for exponents 2/5 (m2) and 1/5 (P) for 
a = 1.5, and 2/7 (1712) and 1/7 (P) for a = 2.5. Insets: the 
compactness index as a function of time in linear-log plots 
for E = 0.35 (a = 1.5) and E = 0.45 (cr = 2.5). 



FIG. 2: Time evolution of the normalized energy distribution 
Ei/E of wave packets in the neighborhood of the initially 
excited site Iq = 500 for o = 2.5 and total energy E — 0.45 
(weak chaos regime, left plot), E = 2.0 (selftrapped regime, 
right plot). Darker regions correspond to higher intensities. 
The disorder realization is the same as in Fig. O 



times, if the initial energy is fixed, and a lowered. For 
low enough initial energies the regime of strong chaos 
can be again avoided. However low initial energies im- 
ply large time scales which are needed to detect any type 
of spreading [9, 11]. Since for a = the equations be- 
come linear again, Anderson localization is restored for 
all times. Therefore, we expect that for fixed initial en- 
ergy and (7 — > 0, the characteristic time scales diverge as 
well. 

Representative examples in the selftrapping and weak 
chaos regimes for a = 1.5 and a = 2.5, are shown in 
Fig. [3] In the regime of weak chaos wave packets ini- 
tially evolve as in the linear case: i. e. they show Ander- 



son localization up to some time Td and both TO2 and P 
remain constant. For t > the wave packet starts to 
grow with TO2 ^ t", P ^ t"/^ (blue curves). The values 
a = 2/5 for a — 1.5 and a = 2/7 for a = 2.5, obtained by 
a theoretical prediction for the weak chaos regime given 
in d, [llj (see also Sect. IIV Al below) . describe quite well 
the numerical data. Increasing the energy shortens the 
time Td (green curves). The evolution of the compactness 
index C for these cases is shown in the insets of Fig. [S] 
For both values of a the compactness index eventually 
oscillates around some constant non-zero value, imply- 
ing that the wave packet does not become more sparse 
when it spreads. In the regime of selftrapping a large 
part of the wave packet remains localized, and therefore 
P is practically constant, while the rest spreads and the 
second moment increases as m2 ^ t" (red curves). The 
time evolution of the energy fraction Ey contained in a 
localization volume V = 20 around the initially excited 
site for a = 2.5 in the cases of immediate subdifFusion 
and of selftrapping is shown in the insets of Fig. [1] 
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IV. SPREADING OF WAVE PACKETS 
A. Theoretical predictions 

In the regime of weak chaos, only a fraction of modes 
in the packet resona ntly interact and facilitate random- 
ization of phases d, [ll[ . The second moment of the 
wave packet increases in time and follows a power law 
of the form m2 ~ t". In these cases, the participation 
number follows the law P ~ t"/^. The dependence of 
the exponent a on the order a of the nonlinearity was 
predicted to be 0, 



The derivation is based on the result that the probability 
V of resonance for a given packet mode depends on the 
energy density e in the packet. For a„f"l'^ > d (strong 
chaos) it follows V ~ 1, while for a^e'^l'^ <^ d (weak 

chaos) we get V sa """^^^ [Tsj . The diffusion rate is 
conjectured to be ~ e'^('P(e))^, which in the case of 
weak chaos, leads to ([TU|) together with m2 ^ 1/e. 

In the regime of strong chaos weget V — 1 and the 
exponent a of m2 ~ is given by [gllisj 



Since the energy density in the packet is decreasing with 
increasing time, the condition for strong chaos will be 
eventually violated, and the spreading will cross over into 
the regime of weak chaos [1^. The crossover duration 
can be very large, and complicates the fitting analysis of 
numerical data. 

If the phases of normal modes are randomized by an 
explicit routine during the evolution of the system, then 
the regime of strong chaos is enforced irrespectively on 
whether the energy density satisfies the criterion of strong 
or weak chaos. In that case the spreading should again 
follow the law pT|) but now for all times [9] . 

The above predictions are expected to hold if a coher- 
ent transfer of energy can be neglected, and only inco- 
herent diffusive transfer is relevant. 

In previous studies with a = 2 and single site excita- 
tions, spreading should start in the regime of weak chaos. 
Corresponding fits of the exponent yield a = 0.33 ± 0.05 
for the KG model [1, in agreement with Eq. ([T0| . 
Mulansky computed spreading exponents for the gDNLS 
model with single site excitations and a = 1,2,4,6 12]. 
Since for a = 2,4,6 strong chaos is avoided, the fitting 
of the dependence m2{t) with a single power law is rea- 
sonable. The corresponding fitted exponents 0.31 ± 0.04 
(ct = 2), 0.18±0.04 (cr = 4) and 0.14±0.05 (cr = 6) agree 
well with the predicted weak chaos result 1/3,1/5,1/7 
from (fTO)) . For a = 1 the initial condition has been 
launched in the regime of strong chaos. A single power 
law fit will therefore not be reasonable. Since the out- 
come is a mixture of first strong and later possibly weak 



chaos, the fitted exponent should be a number which is 
located between the two theoretical values 1/2 and 2/3. 
Indeed, Mulansky reports a number 0.56 ± 0.04. Veksler 
et al. [13 considered short time evolutions of single site 
excitations (up to t ^ 10"^) for gDNLS models. While 
the time window may happen to be too short for con- 
clusive results, the observed increase of fitted exponents 
with increasing /3 for < 2 is possibly also influenced by 
the transition from weak to strong chaos. 



B. Numerical results 

In our study we provide numerical evidence for the 
validity of both predictions ([T0| and ([TT|) , by performing 
extensive simulations of the gKG model for various values 
of (T and for energies E away from the selftrapping regime. 
We consider not only integer but also noninteger values 
of a. 

The used energies E are not too small in order to avoid 
long delays on the onset of spreading, and their values 
cross the boundary between the weak and strong chaos 
regimes as a decreases (Fig.[T|). The transition from weak 
to strong chaos is not an abrupt one. One should keep in 
mind that the curves plotted in Fig. [T] should be consid- 
ered as rough indicators of the borders between different 
regimes, since they are influenced by the characteristics 
of each particular disorder realization. 



1. Computational techniques 

For several values of the nonlinearity order a we com- 
pute the evolution of single site excitations up to a large 
final time tfin for an energy value which excludes self- 
trapping, and for 20 different disorder realizations. For 
each case we verify that the time evolution of m2 and P 
indicate that the wave packet is not selftrapped. Al- 
though we also compute P we chose to analyze and 
present results only for m2, because it grows faster than 
P allowing a more accurate determination of the expo- 
nent a. Typically tfin ranges from tfin — 10^ for small 
cr to t fin — 10^ for larger a for which slower spreading is 
observed. For each realization the m2 (t) is fitted by a t" 
law, the value of the exponent a is determined and the 
average value (a) over the 20 realizations is computed. 
In practice, we perform a linear fit of the log]^QTO2(i), 
logj^Q t values, instead of a nonlinear fit of the actual 
7712 (t), t values. The data used for the fitting lie in the 
time window [logj^Q tini, log^o tfin] whose lower end varies 
from logio t^m = 1 up to log^p = logio " 2, in or- 
der to guarantee that even the smallest window contains 
enough data points (for at least two orders of magnitude 
of t) allowing a reliable evaluation of a. 

If the time evolution of 7712 (t) is well approximated by 
a t" law the averaged value (a) should not depend on the 
initial time tini of the time window. In other words, if 
(a) eventually becomes a constant function of logj^Q tini 



for large enough values of t^i then we consider that t^"^ 
satisfactorily describes the asymptotic (t — oo) behavior 
of m2{t). If on the other hand, (a) does not tend to be- 
come constant as Uni increases, the available numerical 
data cannot be represented reliably by a fit and no 
exponent can be determined. Such an example is seen 
in Fig. Slja) where we see that for a = 0.8 (a) increases 
monotonically as tini increases. The horizontal dashed 
line denotes the theoretically predicted weak chaos ex- 
ponent from Eq. ([TUl) a = 5/9 = 0.556. At the largest 
values of ti^i the exponent is close to the strong chaos 
result a = 5/7 = 0.714 from Eq. ((TTl) (horizontal dotted 
line). 

For a = 2.5 (Fig. IHb)) the values of (a), after some 
initial transient time, seem to saturate to a constant value 
implying that a t" law well approximates the evolution of 
TO2- Actually, (a) eventually attains values close to the 
theoretically predicted weak chaos exponent a — 2/7 = 
0.286, denoted by a horizontal line in Fig.UJ^b). 

A significant dynamical quantity is the minimum 
time T* after which (a) becomes practically constant, 
since it characterizes the onset of the validity of the 
asymptotic approximation m2 ~ t". The value of r* 
is estimated from the numerically evaluated function 
(a) = /(logj^g Uni) as the minimum time such that for 
all Uni > T* the local (numerically estimated) rate 
d{{a)) / d{\ogiQ tini) guarantees a less than 5% change of 
(a) at tfin with respect to its current value at r*. This 
definition has a degree of arbitrariness but it manages 
to capture for all our simulations the time after which 
a fit seems to approximate quite accurately the val- 
ues of 7712 (t). For example for a = 2.5 (Fig. mjb)), we 
find log^o T* = 4, and the corresponding exponent value 
(a) = 0.288 ± 0.058 is very close to the theoretically pre- 
dicted value a = 0.286 from Eq. ([T0|. 

The cases of a = 2.5 and a = 0.8 shown in Fig. 2] are 
two typical examples where an exponent a for TO2 {t) ^ t" 
can or cannot be defined respectively. The time evolution 
of TO2 for three particular disorder realizations for cr = 2.5 
and a = 0.8 is plotted in Fig. \M,c). For a = 2.5 the 
three curves tend to increase according to the theoretical 
prediction given in Eq. ((TO)) (solid line). On the other 
hand, for a = 0.8 m2 seems to increase with an increasing 
rate, not showing a tendency to approach any power law, 
at least up to the final integration time tfin — 10*. In 
this case, TO2 deviates from the weak chaos theoretical 
power law t^^^ prediction given in Eq. ([TU| (dashed line) . 
Although a reliable numerical estimation of a constant 
(a) was not possible from our simulations, the values of 
TO2 for large values of t seem to be close to the strong 
chaos prediction m2 ^ t^^"^ given from Eq. ([TT|) (dotted 
line). The deviation from the power law predictions for 
CT = 0.8 can also be clearly seen by plotting the mean 
value (logj^o "^2) of logj^g ^2{t) over the 20 realizations as 
a function of time (Fig. HJ^d)). Similar results for a — 2.5 
show that the numerical results are in good agreement 
with the weak chaos prediction 1112 ~ t^^^. 
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FIG. 4: (Color online) The numerically obtained mean value 
(a) of the exponent a [ra-2 ~ t") over 20 disorder realiza- 
tions as a function of the initial time ti^i of the fitting win- 
dow having as final time tji^ ~ 10* in linear- log plots for (a) 
a = 0.8, E = 0.25 and (b) a = 2.5, E = 0.45. (c) logn, mz 
versus logj^g t for three difi^erent disorder realizations (denoted 
by (b) blue, (g) green, and (r) red) for cr = 0.8 (dashed up- 
per curves) and a — 2.5 (solid lower curves), (d) The mean 
value (log 7712) of logjQ m2 (with error bars) over 20 disor- 
der realizations versus logj^Q t for a = 0.8 (upper curve) and 
a = 2.5 (lower curve). In all panels straight lines correspond 
to theoretical predictions for the values of a. For a — 0.8 the 
weak chaos exponent a = 5/9 from Eq. (|10|) (dashed lines) 
and the strong chaos exponent a = 5/7 from Eq. (|11|) (dotted 
lines) are plotted, while for a — 2.5 the weak chaos exponent 
a = 2/7 (solid lines) is shown. 



2. Subdiffusive spreading 

We performed extensive simulations for 25 different 
values of the nonlinearity order a in the interval < 
cr < 4. For each a, we followed the evolution of single 
site excitations for 20 different disorder realizations by 
considering an energy value away from the selftrapping 
regime, which allows the immediate subdiffusion of the 
wave packet. For each case we tried to determine the 
exponent a {m,2 ^ t") and the time r* by the above- 
described procedure. Our criterion for determining a 
and T* was satisfied for a — 0.05 and 1.25 < cr < 4. 
For all other tested values of cr (0.1 < cr < 1) a reliable 
value for a was not obtained, since 7712 exhibited behav- 
iors similar to the one seen for ct = 0.8 (Fig. |4]). It is 
possible that for these cases integration up to larger final 
times might allow TO2 to attain its asymptotic power law 
behavior, and permit the estimation of a, but the ex- 
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FIG. 5: (Color online) Exponent a (m2 ^ t") versus the 
nonlinearity order o for plain integration without dephasing 
(filled squares) and for integration with dephasing of NMs 
(filled triangles). Results without dephasing obtained in 
are plotted with empty circle symbols. The theoretically pre- 
dicted functions a — l/{l + a) (weak chaos) and a — 2/{2 + a) 
(strong chaos) are plotted by dashed and solid lines respec- 
tively. Inset: The logarithm of the minimum time r* for which 
the evolution of m2 can be numerically fitted by a function 
of the form t°' versus a for integration with (filled triangles) 
and without (filled squares) dephasing. 



tremely long CPU times needed for such simulations do 
not make them easily feasible. As we can see from Fig.[l] 
all these cases belong to the strong chaos regime. The 
fact that exactly for cr < 2 the intermediate (and possibly 
rather long lasting) regime of strong chaos may set in, is 
a possible explanation for the observed difficulties. 

The computed exponents a for different values of a 
are plotted in Fig. [5] (filled squares) . The values of a ob- 
tained in [13] for gDNLS are also plotted (empty circles). 
These values are very close to our results for a = 2 and 
cr = 4, while for a = 1 a reliable estimate of a was not 
obtained from our simulations. In Fig. [5] the theoreti- 
cally predicted law pUj) is plotted by a dashed line and 
all computed exponents a are in good agreement with 
it. The exponents for a = 1.75, a — 1.5 and a = 1.25 
slightly deviate from this law, possibly due to the infiu- 
ence of the strong chaos regime which exists for < 2. 
The time t* after which m2 is well approximated by 
(inset of Fig. [5|), increases as a approaches zero, imply- 
ing that integration for longer times might be needed in 
order to estimate a for 0.1 < ct < 1. 

From Fig. [T] we expect that for small values of a the 
range of energies E for which we observe immediate sub- 
diffusive spreading in the regime of strong chaos should 
increase. As an example we consider the case of cr = 0.05, 
for which the exponent a was obtained for E — 0.3 
(Fig. [5]). In Fig. [HJa) the evolution of (log^Q (av- 
erage value over 20 realizations) is plotted for E — 0.3 
(green curve), E = 5000.0 (red curve) and E = 50000.0 
(blue curve). For all energies {\0giQm2) increases in a 




FIG. 6: (Color online) (a) mean value (logj^g m2) of logj^Q m2 
over 20 realizations versus logjg t for a = 0.05 and en- 
ergy E = 0.3,5000,50000 [(g) green; (r) red; (b) blue]. 
Straight line guide the eye for the theoretically predicted 
strong chaos exponent a = 0.976. (b) logj^Q 7712 versus 
logjQ t for the same disorder realization with E = 0.03 and 
cr = 0.00001,0.005,0.01,0.05 [(bl) black; (r) red; (g) green; 
(b) blue]. 



similar way following a power law, after some transient 
initial time interval, which is well approximated by the 
theoretically predicted function pT|) t^-^"^^ (dashed line). 
It is evident that the exponent a, which for E = 0.3 was 
estimated to be = 0.93 ± 0.11 at logj^Q r* = 5.4, does 
not depend on the energy value. Note that the error bar 
is too large here to discriminate between the strong chaos 
a — 0.976 and weak chaos a — 0.952 predictions. 

For very small values of cr the dynamics should ap- 
proach the behavior of the linear system (cr = 0), i. e. the 
wave packet should be localized. This tendency is seen 
in Fig. [6l^b) where \ogiQm2{t) is plotted for a disorder 
realization with E = 0.3 in the cases of cr = 0.05, 0.01, 
0.005, 0.00001 (curves from top to bottom). It is evident 
that the time t* increases as cr —>■ 0. In particular, for 
cr = 0.00001 this time was not reached until the end of 
the simulation at t = 10*. 

The order of nonlinearity a infiuences not only the 
spreading rate of wave packets, but also the morphology 
of their profiles. In Fig. [7] we plot the normalized energy 
distributions of initial single site excitations, for different 
cr values in NM (upper plot) and real (lower plot) space. 
Starting from the outer, most extended wave packet we 
plot distributions for cr = 0.05 (black curves), cr = 0.2 
(magenta curves), a ~ 0.8 (red curves), cr = 1.25 (blue 
curves), cr = 2 (green curves) and cr = 3 (brown curves). 
All wave packets were considered for the same disorder 
realization but at different times of their evolution when 
they have the same value of second moment 7712 ~ 10^. 
These times are t = 3.6 x 10^ for cr = 0.05, t = 1.3 x 10^ 
for cr = 0.2, t = 2.5 X 10^ for cr = 0.8, t ^ 1.4 x 10^ for 
cr 1.25, t = 3 X 10'^ for cr = 2 and t = 10^ for cr = 3 and 
increase for a > 0.2 since the spreading becomes slower 
for larger cr. As we have seen in Fig. [51 when cr wave 
packets remain localized for very large time intervals be- 
fore they start to spread. This is why for cr = 0.05 the 
second moment becomes 1712 ~ 10'^ at a larger time than 
in cases with cr = 0.2 and cr = 0.8. 
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FIG. 7: (Color online) Normalized energy distributions in 
NM (upper plot) and real (lower plot) space for a — 
0.05,0.2,0.8,1.25,2.0,3.0 [(bl) black; (m) magenta; (r) red; 
(b) blue; (g) green; (br) brown] at times t = 3.6 x 10^, 1.3 x 
10^ 2.5 X 10^ 1.4 X 10^ 3 x lO^ 10" respectively. The second 
moment of each distribution is m2 ~ 10^ In the upper plot 
the distributions for a — 1.25,2.0 are not clearly visible as 
they are overlapped by the distribution for a — 3.0. 



From the results of Fig.[7]we see that for large enough 
values of a (0.8 < c < 3), the distributions on a loga- 
rithmic scale have a chapeau-like shape consisting of a 
highly excited central part and exponential tails having 
practically the same slope. Contrarily, the distributions 
for (7 = 0.2 and a = 0.05 become more extended having 
different slopes in the tails. 

A characteristic of the NM space distributions in FiglT] 
for a > 0.8 is that they exhibit very large value fluctua- 
tions (up to 10-15 orders of magnitude) in their tails, con- 
trarily to the corresponding distributions in real space. 
Tail NMs are driven by the core of the wave packet, 
and may also interact with neighboring tail NMs. The 
presence of large tail amplitude fluctuations signals that 
neighboring tail NMs do not interact significantly (other- 
wise we would expect a tendency towards equipartition). 
Tail NMs are then excited only by the core; the further 
away they are, the weaker the excitation. But within a 
small tail volume, NMs with larger localization length 
will be more strongly excited than those with smaller lo- 
calization length, hence the large observed fluctuations, 
which on a logarithmic scale are of the order of the rel- 
ative variation of the localization length. Therefore An- 
derson localization is preserved in the tails of the distri- 
butions over very long times (essentially until the given 
tail volume becomes a part of the core). But the NM 
space distributions for a = 0.05 and a — 0.2 exhibit less 
fluctuations in their tail values with respect to the other 
distributions in the upper plot of Fig. [71 implying that 
tail NMs are now interacting with each other on com- 
paratively short time scales and reach a visible level of 



local equipartition. Therefore we observe for these cases 
a destruction of Anderson localization even in the tails 
of the spreading wave packets. 

3. Dephasing 

The assumption that all NMs in a wave packet are 
chaotic leads to a power law increase of m2 whose expo- 
nent a is given by Eq. (|lip . In order to check the validity 
of this prediction we enhanced the wave packet chaotic- 
ity by a periodic dephasing of its NMs. Every 100 time 
units on average 50% of the NMs were randomly chosen 
and the signs of their momenta were changed. In this 
way a faster spreading of the wave packet, with respect 
to its normal evolution, was achieved. We also note that 
a has a similar effect on the shape of wave packets as in 
the case of normal evolution without dephasing, because 
energy distributions in NM and in real space for the same 
7712 values and different cr values have similar profiles to 
the ones shown in Fig. [T] 

Performing a similar numerical analysis (as in the case 
of normal wave packet evolution) we computed the expo- 
nent a of 7)7,2 ~ and the time r* for several values of 
the nonlinearity order a. The obtained values are plotted 
in Fig. [5] by filled triangles. The numerically computed 
exponents are in good agreement with the theoretical pre- 
diction of Eq. ((TT|l (solid line in Fig. [5]). In the case of 
normal wave packet evolution not all NMs in the packet 
are chaotic, because the exponents presented by filled 
squares in Fig.[5]are always smaller that the prediction of 
Eq. pT|) . The dephasing procedure increases drastically 
the chaotic nature of the dynamics since the correspond- 
ing exponents a are quite close, but somewhat smaller, 
than the predicted values given by Eq. pT|) (strong chaos 
regime). 

From the results of Fig. [5] we see that exponents a 
were determined for a larger value interval of a (0.2 < 
(T < 4) with respect to the normal evolution case. Only 
for a = 0.05 and a = 0.1 we were not able to estimate an 
exponent a from the performed numerical simulations. 

In addition, time r* (inset of Fig. O lias always larger 
values with respect to the ones obtained for the nor- 
mal evolution of wave packets, especially for large val- 
ues of a. This means that although dephasing increases 
the chaoticity of the wave packet, a considerably large 
amount of time is needed for the evolution to be charac- 
terized by 7772 t" . 



V. SUMMARY AND DISCUSSION 

We performed extensive numerical simulations of the 
evolution of single site excitations in the gKG model ([3]) 
for different values of the nonlinearity order a. According 
to the analytical treatment presented in [l^ , in this case 
a wave packet could either: a) be selftrapped for large 
enough values of the nonlinearity, i. e. the total energy 
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E of the gKG system, or b) spread subdiffusively for 
smaller values of E. Particularly for energy values not in 
the selftrapping regime, the single site excitation belongs 
either to the weak or the strong chaos regime 

In the weak chaos regime the wave packet spreads sub- 
diffusively and its second moment m2 grows as m2 ~ t". 
The expected dependence of a on cr in this case is given 
by Eq. ([T0|) . The detrapping time tj, after which spread- 
ing stars, increases as E decreases because the system is 
closer to a linear model where no spreading is observed 
due to Anderson localization. In order to be able to ob- 
serve spreading for large time intervals we avoided very 
small energy values. 

If the wave packet is launched in the strong chaos 
regime its subdiffusive spreading is initially characterized 
by an exponent a given from Eq. which is expected 
to eventually cross over to its asymptotic (t — >■ oo) value 
given by Eq. (jlOp . The time at which this crossover starts 
(as well as its duration) could become very large, limiting 
our ability to observe it numerically. 

According to the estimations of fis'l, if the single site 
excitation is not selftrapped then it belongs to the weak 
chaos regime for ct > 2. For a < 2 the existence of the 
strong chaos regime is possible, but not for very small 
energy values where the dynamics is again in the weak 
chaos regime (see Fig. [T|). 

We performed numerical simulations for various val- 
ues of a and for E belonging both to the weak and the 
strong chaos regimes, having care to avoid very small en- 
ergy values where we might encounter very large detrap- 
ping times (Fig. [T]). The curves in Fig. [T]do not define 
exactly the separation between different regimes, instead 
they should be considered as indications of the location 
of border regions between them. The energy values used 
here cross the boundary between the weak and strong 
chaos regimes around the interval 1 ^ cr < 2. 

Our results verify the validity of Eq. ([TU)) for the weak 
chaos regime. In particular, the numerically computed 
exponents a are in very good agreement with the theo- 
retical prediction of Eq. ([TU|) for u > 2, while they exhibit 
an increasing deviation from this prediction for a = 1.75, 
a — 1.5 and a ~ 1.25 respectively (Fig. [5]). For a < I 
m2 grows faster than the weak chaos estimation ([TU]). 
but up to the final times that we performed computa- 
tions we were not able to reliably fit its evolution with 
a power law and compute an exponent a. Nevertheless, 
this problem, as well as the deviation of the computed ex- 
ponents from the theoretical weak chaos prediction ^TU\\ 
for 1.25 < a < 1.75, clearly indicate that in our simula- 
tions the boundary between the weak and strong chaos 
regimes lies in the interval 1 ^ c < 2, being in agreement 
with the theoretical predictions of Fig. [T] 

Regarding the simulations for ct < 1 we note that 
maybe longer integrations could allow for a reliable es- 
timation of a. This is a very hard computational task 
as it requires large CPU times. Since simulations with ct 
values around ct = 1 are located close to the borders of 
different regimes, it is possible that for different disorder 



realizations the wave packet's evolution is a mixture of 
different dynamical behaviors, not allowing the statisti- 
cal analysis to clearly determine a. In addition, for small 
values of ct we observe differences in the wave packet dy- 
namics, which could affect the determination of a, since 
Anderson localization is destroyed also in the tails of the 
wave packets (Fig [7]). 

In order to enforce the wave packet to evolve continu- 
ously in the strong chaos regime we enhanced its chaotic- 
ity by repeatedly performing a dephasing of its NMs. 
When dephasing was applied, the computed exponents a 
were in good agreement with the strong chaos theoreti- 
cal prediction of Eq. pT|) for almost all tested values of 
CT(Fig.[5|. 

Predictions (HO]) and ([TT]) were derived in Refs. [|, [ll] 
for the gDNLS model ([2]) and for integer values of ct. 
Thus, as a final remark we note that our results also 
establish the generality of these predictions since they 
proved to be valid for a different dynamical system (the 
gKG model ([3])) and for noninteger values of ct. 
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Appendix A: Frequency shift of nonlinear oscillators 

Let us consider a nonlinear oscillator described by the 
Hamiltonian function 



9 9 2 

2 2 



/3 



CT + 2' 



,|<T+2 



(Al) 



with CT > 0. For a given value of the energy E > the 
oscillator's period T is 



T = 4 



dx 



2E - LO^X^ - 



(A2) 



2^2 



where x is the positive root of equation uj x 
2E. The change of variable x 

^2Evy[vi) /uj, with 



2/3 
cr+2 I 



2/3 {2Ey/^ 



(A3) 



CT + 2 

and y{ri) being the positive number satisfying 

y{vf + v\y{vW^^^h (A4) 

transforms Eq. (jA2p to 

4y(??) dv 



Tiv) 



U! Jq ^\ — y(ry)2-y2 _ jjy(^j^'^cr+2yCr+2 



(A5) 
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For /? = we get t] = 0, y = 1 and Eq. (jAip cor- 
responds to an harmonic oscillator with period T(0) = 
2tt/uj and frequency 17(0) = uj. In order to estimate the 
change of the squared frequency Sft^ = ^^(77) — il^{0), 
we note that the period of the oscillator for 77 > is 
T{ri) « T{0) + rjT'iQ), at a first order approximation in 
?7, with ' ' ' denoting derivative with respect to r/. Dif- 
ferentiating Eq. (IA4[) we find y'(0) — —1/2. Using this 
result, and differentiating Eq. (jASp according to the Leib- 
niz integral rule we get 



r'(o) = - 



27r 



1 



7r7o (1-^2)3/2 
Then, frequency il(?7) is given by 



dv 



(A6) 



Since Q 



(1-^2)3/2 

frequency shift is 



r((^+3)/2) 
r((a+2)/2)' 



tT+4 



r(^ 



V^(a + 2)r(^) 



the squared 



(A8) 



27r 



1 w2-t,-+2 



^^''(2^7r7o (l-?;2)3/2 



dv 



(A7) Eq. © is derived from Eq. ([M]) for /3 = 1 and = 
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